Nesta lista, vamos trabalhar com o seguinte processo estocástico AR(1):

\[ z_t = \rho z_{t-1} + \varepsilon_t , \]

com \(\varepsilon_t \sim N(\mu,\sigma^2)\). Inicialmente, assumiremos que \(\mu=0\), \(\rho = 0.95\) e \(\sigma = 0.007\), calibração de Cooley & Prescott (1995).


Questão 1

Discretize o processo acima usando o método de Tauchen (1986). Use 9 pontos.


Ao discretizar um processo que assume valores contínuos, é necessário determinar o intervalo do grid, o número de pontos, e sua localização. Além disso, é necessário determinar probabiliades condicionais de transição entre os estados de forma consistente com o processo original.

Pelo método de Tauchen, o intervalo de um grid de \(N\) pontos é determinado por \([\mu - \theta_N, \mu + \theta_N]\), em que \(\theta_N = m \frac{\sigma}{\sqrt{1-\rho^2}}\) é o desvio padrão do processo estocástico, multiplicado por um parâmetro de escala \(m\). Os demais \(N-2\) pontos são distribuídos uniformemente no intervalo.

# Parâmetros Iniciais
N     = 9
rho   = 0.95
sigma = 0.007
mu    = 0
m     = 3
tauchen_grid = function(N, m, sigma, rho) {
  thetaN = m * sigma / sqrt(1-rho^2)  
  seq(-thetaN, thetaN, length.out = N)
}
thetaT = tauchen_grid(N, m, sigma, rho)
as.matrix(thetaT)
             [,1]
 [1,] -0.06725382
 [2,] -0.05044037
 [3,] -0.03362691
 [4,] -0.01681346
 [5,]  0.00000000
 [6,]  0.01681346
 [7,]  0.03362691
 [8,]  0.05044037
 [9,]  0.06725382

Seja \(\theta_{i,t}\) um estado no período \(t\), e \(\theta_{j, t+1}\) um estado no período \(t+1\). No método de Tauchen, as probabilidades de transição dos pontos interiores do grid são dadas por

\[ p_{i,j} = \Phi \left( \frac{\theta_j + \Delta \theta/2 - (1-\rho) \mu - \rho \theta_i}{\sigma} \right) - \Phi \left(\frac{\theta_j - \Delta \theta/2 - (1-\rho) \mu - \rho \theta_i}{\sigma}\right)\]

em que \(\Phi\) é a CDF da \(N(\mu, \sigma^2)\). Nos cantos do grid, temos que

\[ p_{i,1} = \Phi \left(\frac{\theta_1 - (1-\rho) \mu - \rho \theta_i + \Delta \theta/2}{\sigma} \right) \]

\[ p_{i,N} = 1 - \Phi \left(\frac{\theta_N - (1-\rho) \mu - \rho \theta_i - \Delta \theta/2}{\sigma} \right).\]

tauchen_P = function(grid, rho, sigma, mu) {
  N = length(grid)
  delta = (max(grid) - min(grid)) / (N-1)
  PT = matrix(NA_real_, N, N)
  for(i in 1:N) {
    for(j in 1:N) {
      if(j == 1) {
        PT[i,j] = pnorm( (grid[1] - (1-rho)*mu - rho*grid[i] + delta/2) / sigma )
      } else if(j==N) {
        PT[i,j] = 1 - pnorm( (grid[N] - (1-rho)*mu - rho*grid[i] - delta/2) / sigma )
      } else {
        PT[i,j] = 
          pnorm((grid[j] + delta/2 - (1-rho)*mu - rho*grid[i]) / sigma) -
          pnorm((grid[j] - delta/2 - (1-rho)*mu - rho*grid[i]) / sigma)
      }
    }
  }
  return(PT)
}
PT = tauchen_P(thetaT, rho, sigma, mu)
round(PT, 4)
        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
 [1,] 0.7644 0.2347 0.0009 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000
 [2,] 0.0592 0.7405 0.1997 0.0006 0.0000 0.0000 0.0000 0.0000 0.0000
 [3,] 0.0001 0.0747 0.7569 0.1679 0.0004 0.0000 0.0000 0.0000 0.0000
 [4,] 0.0000 0.0001 0.0931 0.7669 0.1396 0.0002 0.0000 0.0000 0.0000
 [5,] 0.0000 0.0000 0.0002 0.1147 0.7702 0.1147 0.0002 0.0000 0.0000
 [6,] 0.0000 0.0000 0.0000 0.0002 0.1396 0.7669 0.0931 0.0001 0.0000
 [7,] 0.0000 0.0000 0.0000 0.0000 0.0004 0.1679 0.7569 0.0747 0.0001
 [8,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0006 0.1997 0.7405 0.0592
 [9,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0009 0.2347 0.7644
rowSums(PT) # Linhas da matriz devem somar 1.
[1] 1 1 1 1 1 1 1 1 1


Questão 2

Discretize o processo acima usando o método de Rouwenhorst. Use 9 pontos.


A definição do grid no método de Rouwenhorst é um caso particular do método de Tauchen, com \(m = \sqrt{N-1}\).

thetaR = tauchen_grid(N, m = sqrt(N-1), sigma, rho)
thetaR
[1] -0.06340751 -0.04755564 -0.03170376 -0.01585188  0.00000000  0.01585188  0.03170376  0.04755564  0.06340751

A matriz de transição \(P\), porém, é obtida recursivamente.

rouwen_P = function(N, rho) {
  p = (1+rho)/2
  P_init = matrix(c(p, 1-p, 1-p, p), 2, 2, byrow = TRUE)
  for(i in 3:N) {
    PR = 
      rbind(cbind(p*P_init, 0), 0) +
      rbind(cbind(0, (1-p)*P_init), 0) +
      rbind(0, cbind((1-p)*P_init, 0)) +
      rbind(0, cbind(0, p*P_init))
    PR = PR / rowSums(PR) # normalização
    P_init = PR
  }
  
  return(PR)
}
PR = rouwen_P(N, rho)
round(PR, 4)
        [,1]   [,2]   [,3]   [,4]   [,5]   [,6]   [,7]   [,8]   [,9]
 [1,] 0.8167 0.1675 0.0150 0.0008 0.0000 0.0000 0.0000 0.0000 0.0000
 [2,] 0.0209 0.8204 0.1469 0.0113 0.0005 0.0000 0.0000 0.0000 0.0000
 [3,] 0.0005 0.0420 0.8231 0.1261 0.0081 0.0003 0.0000 0.0000 0.0000
 [4,] 0.0000 0.0016 0.0630 0.8247 0.1051 0.0054 0.0001 0.0000 0.0000
 [5,] 0.0000 0.0001 0.0032 0.0841 0.8253 0.0841 0.0032 0.0001 0.0000
 [6,] 0.0000 0.0000 0.0001 0.0054 0.1051 0.8247 0.0630 0.0016 0.0000
 [7,] 0.0000 0.0000 0.0000 0.0003 0.0081 0.1261 0.8231 0.0420 0.0005
 [8,] 0.0000 0.0000 0.0000 0.0000 0.0005 0.0113 0.1469 0.8204 0.0209
 [9,] 0.0000 0.0000 0.0000 0.0000 0.0000 0.0008 0.0150 0.1675 0.8167


Questão 03

Simule o processo contínuo para 10000 períodos. Faça o mesmo para os processos discretizados (lembre-se de usar as mesmas realizações para os choques). Compare os caminhos para cada processo (gráficos serão úteis aqui). Se eles não estiverem muito próximos, utilize mais pontos.


set.seed(96452)
time_span = 10000
eps = rnorm(time_span, mu, sigma)
# Processo contínuo
zc = rep(0, time_span)
for(t in 2:time_span)
  zc[t] = rho * zc[t-1] + eps[t]
# Funcão para gerar processos
gen_process = function(grid, P, eps, mu, sigma) {
  time_span = length(eps)
  z = rep(0, time_span)
  for(t in 2:time_span) {
    i = which(grid == z[t-1])
    j = sum(cumsum(P[i, ]) < pnorm(eps[t], mu, sigma)) + 1
    z[t] = grid[j]
  }
  return(z)  
}
zt = gen_process(thetaT, PT, eps, mu, sigma)
zr = gen_process(thetaR, PR, eps, mu, sigma)

Os gráficos a seguir exibem o processo contínuo e as discretizações de Tauchen e de Rouwenhorst para 10.000 períodos e para os 500 primeiros períodos:

Visualmente, a discretização de Tauchen funciona melhor com os parâmetros dados. A discretização de Tauchen captura melhor a amplitude do processo sem perder precisão. Ambos apresentam, contudo um RMSE semelhante em relação ao processo contínuo:

RMSE = function(x,y) sqrt(mean((x-y)^2))
cat("RMSE:\n")
cbind("Tauchen" = RMSE(zc,zt), "Rouwenhorst" = RMSE(zc,zr))

Questão 4

Estime processos AR(1) com base nos dados simulados, tanto a partir do Tauchen quanto o de Rouwenhorst. Quão próximo eles estão do processo gerador de dados real? Se eles não estiverem muito próximos, utilize mais pontos.

fit_zt = lm(zt ~ shift(zt) + 0)
fit_zr = lm(zr ~ shift(zr) + 0)
coefs = c(fit_zt$coefficients,
          fit_zr$coefficients)
sigmas = c(summary(fit_zt)$sigma,
           summary(fit_zr)$sigma)
tbl = data.frame(method = c("Tauchen 9 pts", 
                            "Rouwenhorst 9 pts"), 
                 rho = coefs, sigma = sigmas)
tbl

Pelo método de Tauchen com 9 pontos no grid, estimamos \(\hat{\sigma} = 0.008\), valor relativamente distante do \(\sigma\) verdadeiro. Aumentando o número de pontos no grid para 51, obtemos uma estimativa semelhante à obtida para o método de Rouwenhorst com 9 pontos.

thetaT = tauchen_grid(N=51, m=3, sigma, rho)
PT = tauchen_P(thetaT, rho, sigma, mu)
ztl = gen_process(thetaT, PT, eps, mu, sigma)
fit_ztl = lm(ztl ~ shift(ztl) + 0)
rbind(tbl, cbind(method = "Tauchen 25 pts", rho = fit_ztl$coefficients, sigma = summary(fit_ztl)$sigma))

Questão 5

Vamos refazer os exercícios acima com \(\rho = 0.99\).

O método de Rouwenhorst funciona melhor para \(\rho = 0.99\) com 9 pontos no grid. Obtemos os seguintes valores de \(\hat{\rho}\) e \(\hat{\sigma}\):

LS0tDQp0aXRsZTogIk3DqXRvZG9zIE51bcOpcmljb3MgMjAxODxicj5MaXN0YSAwMSINCm91dHB1dDoNCiAgaHRtbF9kb2N1bWVudDoNCiAgICBkZl9wcmludDogcGFnZWQNCiAgaHRtbF9ub3RlYm9vazoNCiAgICB0aGVtZTogZmxhdGx5DQogIHBkZl9kb2N1bWVudDogZGVmYXVsdA0KZWRpdG9yX29wdGlvbnM6DQogIGNodW5rX291dHB1dF90eXBlOiBpbmxpbmUNCi0tLQ0KDQpgYGB7ciBpbmNsdWRlPUZBTFNFfQ0KbGlicmFyeShkYXRhLnRhYmxlKQ0KbGlicmFyeSh0aWN0b2MpDQpsaWJyYXJ5KGdncGxvdDIpDQpsaWJyYXJ5KHBsb3RseSkNCm9wdGlvbnMoc2NpcGVuID0gMTApDQpgYGANCg0KPGJyPjxicj4NCg0KTmVzdGEgbGlzdGEsIHZhbW9zIHRyYWJhbGhhciBjb20gbyBzZWd1aW50ZSBwcm9jZXNzbyBlc3RvY8Ohc3RpY28gQVIoMSk6DQoNCiQkIHpfdCA9IFxyaG8gel97dC0xfSArIFx2YXJlcHNpbG9uX3QgLCAkJA0KDQpjb20gJFx2YXJlcHNpbG9uX3QgXHNpbSBOKFxtdSxcc2lnbWFeMikkLiBJbmljaWFsbWVudGUsIGFzc3VtaXJlbW9zIHF1ZSANCiRcbXU9MCQsICRccmhvID0gMC45NSQgZSAkXHNpZ21hID0gMC4wMDckLCBjYWxpYnJhw6fDo28gZGUgQ29vbGV5ICYgUHJlc2NvdHQgKDE5OTUpLg0KDQo8YnI+DQoNCiMjIyBRdWVzdMOjbyAxDQojIyMjIERpc2NyZXRpemUgbyBwcm9jZXNzbyBhY2ltYSB1c2FuZG8gbyBtw6l0b2RvIGRlIFRhdWNoZW4gKDE5ODYpLiBVc2UgOSBwb250b3MuDQoNCjxicj4NCg0KQW8gZGlzY3JldGl6YXIgdW0gcHJvY2Vzc28gcXVlIGFzc3VtZSB2YWxvcmVzIGNvbnTDrW51b3MsIMOpIG5lY2Vzc8OhcmlvIGRldGVybWluYXINCm8gaW50ZXJ2YWxvIGRvIGdyaWQsIG8gbsO6bWVybyBkZSBwb250b3MsIGUgc3VhIGxvY2FsaXphw6fDo28uIEFsw6ltIGRpc3NvLCDDqQ0KbmVjZXNzw6FyaW8gZGV0ZXJtaW5hciBwcm9iYWJpbGlhZGVzIGNvbmRpY2lvbmFpcyBkZSB0cmFuc2nDp8OjbyBlbnRyZSBvcyANCmVzdGFkb3MgZGUgZm9ybWEgY29uc2lzdGVudGUgY29tIG8gcHJvY2Vzc28gb3JpZ2luYWwuDQoNClBlbG8gbcOpdG9kbyBkZSBUYXVjaGVuLCBvIGludGVydmFsbyBkZSB1bSBncmlkIGRlICROJCBwb250b3Mgw6kgZGV0ZXJtaW5hZG8gcG9yIA0KJFtcbXUgLSBcdGhldGFfTiwgXG11ICsgXHRoZXRhX05dJCwgZW0gcXVlICRcdGhldGFfTiA9IG0gXGZyYWN7XHNpZ21hfXtcc3FydHsxLVxyaG9eMn19JA0Kw6kgbyBkZXN2aW8gcGFkcsOjbyBkbyBwcm9jZXNzbyBlc3RvY8Ohc3RpY28sIG11bHRpcGxpY2FkbyBwb3IgdW0gcGFyw6JtZXRybyBkZSBlc2NhbGEgJG0kLiANCk9zIGRlbWFpcyAkTi0yJCBwb250b3Mgc8OjbyBkaXN0cmlidcOtZG9zIHVuaWZvcm1lbWVudGUgbm8gaW50ZXJ2YWxvLg0KDQpgYGB7cn0NCiMgUGFyw6JtZXRyb3MgSW5pY2lhaXMNCk4gICAgID0gOQ0KcmhvICAgPSAwLjk1DQpzaWdtYSA9IDAuMDA3DQptdSAgICA9IDANCm0gICAgID0gMw0KDQp0YXVjaGVuX2dyaWQgPSBmdW5jdGlvbihOLCBtLCBzaWdtYSwgcmhvKSB7DQogIHRoZXRhTiA9IG0gKiBzaWdtYSAvIHNxcnQoMS1yaG9eMikgIA0KICBzZXEoLXRoZXRhTiwgdGhldGFOLCBsZW5ndGgub3V0ID0gTikNCn0NCg0KdGhldGFUID0gdGF1Y2hlbl9ncmlkKE4sIG0sIHNpZ21hLCByaG8pDQphcy5tYXRyaXgodGhldGFUKQ0KYGBgDQoNCg0KU2VqYSAkXHRoZXRhX3tpLHR9JCB1bSBlc3RhZG8gbm8gcGVyw61vZG8gJHQkLCBlICRcdGhldGFfe2osIHQrMX0kIHVtIGVzdGFkbyBubw0KcGVyw61vZG8gJHQrMSQuIE5vIG3DqXRvZG8gZGUgVGF1Y2hlbiwgYXMgcHJvYmFiaWxpZGFkZXMgZGUgdHJhbnNpw6fDo28gZG9zIHBvbnRvcw0KaW50ZXJpb3JlcyBkbyBncmlkIHPDo28gZGFkYXMgcG9yDQoNCiQkIHBfe2ksan0gPSANCiBcUGhpIFxsZWZ0KCBcZnJhY3tcdGhldGFfaiArIFxEZWx0YSBcdGhldGEvMiAtICgxLVxyaG8pIFxtdSAtIFxyaG8gXHRoZXRhX2l9e1xzaWdtYX0gXHJpZ2h0KQ0KIC0gXFBoaSBcbGVmdChcZnJhY3tcdGhldGFfaiAtIFxEZWx0YSBcdGhldGEvMiAtICgxLVxyaG8pIFxtdSAtIFxyaG8gXHRoZXRhX2l9e1xzaWdtYX1ccmlnaHQpJCQgDQoNCmVtIHF1ZSAkXFBoaSQgw6kgYSBDREYgZGEgJE4oXG11LCBcc2lnbWFeMikkLiBOb3MgY2FudG9zIGRvIGdyaWQsIHRlbW9zIHF1ZQ0KDQokJCBwX3tpLDF9ID0gDQogXFBoaSBcbGVmdChcZnJhY3tcdGhldGFfMSAtICgxLVxyaG8pIFxtdSAtIFxyaG8gXHRoZXRhX2kgKyBcRGVsdGEgXHRoZXRhLzJ9e1xzaWdtYX0gXHJpZ2h0KSAkJA0KDQokJCBwX3tpLE59ID0gMSAtIFxQaGkgXGxlZnQoXGZyYWN7XHRoZXRhX04gLSAoMS1ccmhvKSBcbXUgLSBccmhvIFx0aGV0YV9pIC0gXERlbHRhIFx0aGV0YS8yfXtcc2lnbWF9IFxyaWdodCkuJCQgDQoNCmBgYHtyIGVjaG89VFJVRX0NCnRhdWNoZW5fUCA9IGZ1bmN0aW9uKGdyaWQsIHJobywgc2lnbWEsIG11KSB7DQogIE4gPSBsZW5ndGgoZ3JpZCkNCiAgZGVsdGEgPSAobWF4KGdyaWQpIC0gbWluKGdyaWQpKSAvIChOLTEpDQogIFBUID0gbWF0cml4KE5BX3JlYWxfLCBOLCBOKQ0KICBmb3IoaSBpbiAxOk4pIHsNCiAgICBmb3IoaiBpbiAxOk4pIHsNCiAgICAgIGlmKGogPT0gMSkgew0KICAgICAgICBQVFtpLGpdID0gcG5vcm0oIChncmlkWzFdIC0gKDEtcmhvKSptdSAtIHJobypncmlkW2ldICsgZGVsdGEvMikgLyBzaWdtYSApDQogICAgICB9IGVsc2UgaWYoaj09Tikgew0KICAgICAgICBQVFtpLGpdID0gMSAtIHBub3JtKCAoZ3JpZFtOXSAtICgxLXJobykqbXUgLSByaG8qZ3JpZFtpXSAtIGRlbHRhLzIpIC8gc2lnbWEgKQ0KICAgICAgfSBlbHNlIHsNCiAgICAgICAgUFRbaSxqXSA9IA0KICAgICAgICAgIHBub3JtKChncmlkW2pdICsgZGVsdGEvMiAtICgxLXJobykqbXUgLSByaG8qZ3JpZFtpXSkgLyBzaWdtYSkgLQ0KICAgICAgICAgIHBub3JtKChncmlkW2pdIC0gZGVsdGEvMiAtICgxLXJobykqbXUgLSByaG8qZ3JpZFtpXSkgLyBzaWdtYSkNCiAgICAgIH0NCiAgICB9DQogIH0NCiAgcmV0dXJuKFBUKQ0KfQ0KDQpQVCA9IHRhdWNoZW5fUCh0aGV0YVQsIHJobywgc2lnbWEsIG11KQ0Kcm91bmQoUFQsIDQpDQpyb3dTdW1zKFBUKSAjIExpbmhhcyBkYSBtYXRyaXogZGV2ZW0gc29tYXIgMS4NCmBgYA0KPGJyPg0KDQojIyMgUXVlc3TDo28gMg0KIyMjIyBEaXNjcmV0aXplIG8gcHJvY2Vzc28gYWNpbWEgdXNhbmRvIG8gbcOpdG9kbyBkZSBSb3V3ZW5ob3JzdC4gVXNlIDkgcG9udG9zLg0KDQo8YnI+DQoNCkEgZGVmaW5pw6fDo28gZG8gZ3JpZCBubyBtw6l0b2RvIGRlIFJvdXdlbmhvcnN0IMOpIHVtIGNhc28gcGFydGljdWxhciBkbyBtw6l0b2RvIGRlIA0KVGF1Y2hlbiwgY29tICRtID0gXHNxcnR7Ti0xfSQuIA0KDQpgYGB7cn0NCnRoZXRhUiA9IHRhdWNoZW5fZ3JpZChOLCBtID0gc3FydChOLTEpLCBzaWdtYSwgcmhvKQ0KdGhldGFSDQpgYGANCg0KQSBtYXRyaXogZGUgdHJhbnNpw6fDo28gJFAkLCBwb3LDqW0sIMOpIG9idGlkYSByZWN1cnNpdmFtZW50ZS4gDQoNCmBgYHtyfQ0Kcm91d2VuX1AgPSBmdW5jdGlvbihOLCByaG8pIHsNCiAgcCA9ICgxK3JobykvMg0KICBQX2luaXQgPSBtYXRyaXgoYyhwLCAxLXAsIDEtcCwgcCksIDIsIDIsIGJ5cm93ID0gVFJVRSkNCg0KICBmb3IoaSBpbiAzOk4pIHsNCiAgICBQUiA9IA0KICAgICAgcmJpbmQoY2JpbmQocCpQX2luaXQsIDApLCAwKSArDQogICAgICByYmluZChjYmluZCgwLCAoMS1wKSpQX2luaXQpLCAwKSArDQogICAgICByYmluZCgwLCBjYmluZCgoMS1wKSpQX2luaXQsIDApKSArDQogICAgICByYmluZCgwLCBjYmluZCgwLCBwKlBfaW5pdCkpDQogICAgUFIgPSBQUiAvIHJvd1N1bXMoUFIpICMgbm9ybWFsaXphw6fDo28NCiAgICBQX2luaXQgPSBQUg0KICB9DQogIA0KICByZXR1cm4oUFIpDQp9DQoNClBSID0gcm91d2VuX1AoTiwgcmhvKQ0Kcm91bmQoUFIsIDQpDQpgYGANCjxicj4NCg0KIyMjIFF1ZXN0w6NvIDAzDQojIyMjIFNpbXVsZSBvIHByb2Nlc3NvIGNvbnTDrW51byBwYXJhIDEwMDAwIHBlcsOtb2Rvcy4gRmHDp2EgbyBtZXNtbyBwYXJhIG9zIHByb2Nlc3NvcyBkaXNjcmV0aXphZG9zIChsZW1icmUtc2UgZGUgdXNhciBhcyBtZXNtYXMgcmVhbGl6YcOnw7VlcyBwYXJhIG9zIGNob3F1ZXMpLiBDb21wYXJlIG9zIGNhbWluaG9zIHBhcmEgY2FkYSBwcm9jZXNzbyAoZ3LDoWZpY29zIHNlcsOjbyDDunRlaXMgYXF1aSkuIFNlIGVsZXMgbsOjbyBlc3RpdmVyZW0gbXVpdG8gcHLDs3hpbW9zLCB1dGlsaXplIG1haXMgcG9udG9zLg0KDQo8YnI+DQoNCmBgYHtyfQ0Kc2V0LnNlZWQoOTY0NTIpDQp0aW1lX3NwYW4gPSAxMDAwMA0KZXBzID0gcm5vcm0odGltZV9zcGFuLCBtdSwgc2lnbWEpDQoNCiMgUHJvY2Vzc28gY29udMOtbnVvDQp6YyA9IHJlcCgwLCB0aW1lX3NwYW4pDQpmb3IodCBpbiAyOnRpbWVfc3BhbikNCiAgemNbdF0gPSByaG8gKiB6Y1t0LTFdICsgZXBzW3RdDQoNCiMgRnVuY8OjbyBwYXJhIGdlcmFyIHByb2Nlc3Nvcw0KZ2VuX3Byb2Nlc3MgPSBmdW5jdGlvbihncmlkLCBQLCBlcHMsIG11LCBzaWdtYSkgew0KICB0aW1lX3NwYW4gPSBsZW5ndGgoZXBzKQ0KICB6ID0gcmVwKDAsIHRpbWVfc3BhbikNCiAgZm9yKHQgaW4gMjp0aW1lX3NwYW4pIHsNCiAgICBpID0gd2hpY2goZ3JpZCA9PSB6W3QtMV0pDQogICAgaiA9IHN1bShjdW1zdW0oUFtpLCBdKSA8IHBub3JtKGVwc1t0XSwgbXUsIHNpZ21hKSkgKyAxDQogICAgelt0XSA9IGdyaWRbal0NCiAgfQ0KICByZXR1cm4oeikgIA0KfQ0KDQp6dCA9IGdlbl9wcm9jZXNzKHRoZXRhVCwgUFQsIGVwcywgbXUsIHNpZ21hKQ0KenIgPSBnZW5fcHJvY2Vzcyh0aGV0YVIsIFBSLCBlcHMsIG11LCBzaWdtYSkNCmBgYA0KDQpPcyBncsOhZmljb3MgYSBzZWd1aXIgZXhpYmVtIG8gcHJvY2Vzc28gY29udMOtbnVvIGUgYXMgZGlzY3JldGl6YcOnw7VlcyBkZSBUYXVjaGVuIA0KZSBkZSBSb3V3ZW5ob3JzdCBwYXJhIDEwLjAwMCBwZXLDrW9kb3MgZSBwYXJhIG9zIDUwMCBwcmltZWlyb3MgcGVyw61vZG9zOg0KDQpgYGB7ciBlY2hvPUZBTFNFfQ0KdGJsID0gZGF0YS5mcmFtZSh0ID0gMTp0aW1lX3NwYW4sIHpjLCB6dCwgenIpDQoNCnAxID0gZ2dwbG90KHRibCwgYWVzKHRibCR0KSkgKw0KICBnZW9tX2xpbmUoYWVzKHkgPSB0YmwkemMsIGNvbCA9ICJDb250aW51b3VzIiksIHNpemUgPSAwLjIpICsNCiAgZ2VvbV9saW5lKGFlcyh5ID0gdGJsJHp0LCAgY29sID0gIlRhdWNoZW4gRGlzY3JldGl6YXRpb24iKSwgc2l6ZSA9IDAuNSkgKw0KICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSB0aGV0YVQsIHNpemUgPSAwLjEsIGxpbmV0eXBlID0gMywgYWxwaGEgPSAwLjUpICsNCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcyA9IHJvdW5kKHRoZXRhVCwgNCksIGxpbWl0cyA9IGMoLS4xLCAuMSkpICsNCiAgc2NhbGVfY29sb3JfbWFudWFsKG5hbWUgPSAiUHJvY2VzcyIsIHZhbHVlcyA9IGMoIkNvbnRpbnVvdXMiPSAxLCAiVGF1Y2hlbiBEaXNjcmV0aXphdGlvbiI9IDQpKSArDQogIHRoZW1lX2NsYXNzaWMoKSArDQogIHlsYWIoInoiKSArDQogIHhsYWIoInQiKQ0KDQp0YmwyID0gdGJsWzE6NTAwLCBdDQoNCnAyID0gZ2dwbG90KHRibDIsIGFlcyh0YmwyJHQpKSArDQogIGdlb21fbGluZShhZXMoeSA9IHRibDIkemMsIGNvbCA9ICJDb250aW51b3VzIiksIHNpemUgPSAwLjIpICsNCiAgZ2VvbV9saW5lKGFlcyh5ID0gdGJsMiR6dCwgIGNvbCA9ICJUYXVjaGVuIERpc2NyZXRpemF0aW9uIiksIHNpemUgPSAwLjUpICsNCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gdGhldGFULCBzaXplID0gMC4xLCBsaW5ldHlwZSA9IDMsIGFscGhhID0gMC41KSArDQogIHNjYWxlX3lfY29udGludW91cyhicmVha3MgPSByb3VuZCh0aGV0YVQsIDQpLCBsaW1pdHMgPSBjKC0uMSwgLjEpKSArDQogIHNjYWxlX2NvbG9yX21hbnVhbChuYW1lID0gIlByb2Nlc3MiLCB2YWx1ZXMgPSBjKCJDb250aW51b3VzIj0gMSwgIlRhdWNoZW4gRGlzY3JldGl6YXRpb24iPSA0KSkgKw0KICB0aGVtZV9jbGFzc2ljKCkgKw0KICB5bGFiKCJ6IikgKw0KICB4bGFiKCJ0IikNCg0KcGwxID0gZ2dwbG90bHkocDEsIHRvb2x0aXAgPSBOVUxMKQ0KcGwyID0gZ2dwbG90bHkocDIsIHRvb2x0aXAgPSBOVUxMKQ0KDQpwbCA9IHN1YnBsb3QocGwxLCBwbDIsIG5yb3dzID0gMikNCnBsJGhlaWdodCA9IDYwMA0KcGwkd2lkdGggPSAxLjUgKiBwbCRoZWlnaHQNCg0KcGwNCmBgYA0KDQpgYGB7ciBlY2hvPUZBTFNFfQ0KdGJsID0gZGF0YS5mcmFtZSh0ID0gMTp0aW1lX3NwYW4sIHpjLCB6dCwgenIpDQoNCnAxID0gZ2dwbG90KHRibCwgYWVzKHRibCR0KSkgKw0KICBnZW9tX2xpbmUoYWVzKHkgPSB0YmwkemMsIGNvbCA9ICJDb250aW51b3VzIiksIHNpemUgPSAwLjIpICsNCiAgZ2VvbV9saW5lKGFlcyh5ID0gdGJsJHpyLCAgY29sID0gIlJvdXdlbmhvcnN0IERpc2NyZXRpemF0aW9uIiksIHNpemUgPSAwLjUpICsNCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gdGhldGFULCBzaXplID0gMC4xLCBsaW5ldHlwZSA9IDMsIGFscGhhID0gMC41KSArDQogIHNjYWxlX3lfY29udGludW91cyhicmVha3MgPSByb3VuZCh0aGV0YVQsIDQpLCBsaW1pdHMgPSBjKC0uMSwgLjEpKSArDQogIHNjYWxlX2NvbG9yX21hbnVhbChuYW1lID0gIlByb2Nlc3MiLCB2YWx1ZXMgPSBjKCJDb250aW51b3VzIj0gMSwgIlJvdXdlbmhvcnN0IERpc2NyZXRpemF0aW9uIj0gImJyb3duMyIpKSArDQogIHRoZW1lX2NsYXNzaWMoKSArDQogIHlsYWIoInoiKSArDQogIHhsYWIoInQiKQ0KDQp0YmwyID0gdGJsWzE6NTAwLCBdDQoNCnAyID0gZ2dwbG90KHRibDIsIGFlcyh0YmwyJHQpKSArDQogIGdlb21fbGluZShhZXMoeSA9IHRibDIkemMsIGNvbCA9ICJDb250aW51b3VzIiksIHNpemUgPSAwLjIpICsNCiAgZ2VvbV9saW5lKGFlcyh5ID0gdGJsMiR6ciwgIGNvbCA9ICJSb3V3ZW5ob3JzdCBEaXNjcmV0aXphdGlvbiIpLCBzaXplID0gMC41KSArDQogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IHRoZXRhVCwgc2l6ZSA9IDAuMSwgbGluZXR5cGUgPSAzLCBhbHBoYSA9IDAuNSkgKw0KICBzY2FsZV95X2NvbnRpbnVvdXMoYnJlYWtzID0gcm91bmQodGhldGFULCA0KSwgbGltaXRzID0gYygtLjEsIC4xKSkgKw0KICBzY2FsZV9jb2xvcl9tYW51YWwobmFtZSA9ICJQcm9jZXNzIiwgdmFsdWVzID0gYygiQ29udGludW91cyI9IDEsICJSb3V3ZW5ob3JzdCBEaXNjcmV0aXphdGlvbiI9ICJicm93bjMiKSkgKw0KICB0aGVtZV9jbGFzc2ljKCkgKw0KICB5bGFiKCJ6IikgKw0KICB4bGFiKCJ0IikNCg0KcGwxID0gZ2dwbG90bHkocDEsIHRvb2x0aXAgPSBOVUxMKQ0KcGwyID0gZ2dwbG90bHkocDIsIHRvb2x0aXAgPSBOVUxMKQ0KDQpwbCA9IHN1YnBsb3QocGwxLCBwbDIsIG5yb3dzID0gMikNCnBsJGhlaWdodCA9IDYwMA0KcGwkd2lkdGggPSAxLjUgKiBwbCRoZWlnaHQNCg0KcGwNCmBgYA0KDQpWaXN1YWxtZW50ZSwgYSBkaXNjcmV0aXphw6fDo28gZGUgVGF1Y2hlbiBmdW5jaW9uYSBtZWxob3IgY29tIG9zIHBhcsOibWV0cm9zDQpkYWRvcy4gQSBkaXNjcmV0aXphw6fDo28gZGUgVGF1Y2hlbiBjYXB0dXJhIG1lbGhvciBhIGFtcGxpdHVkZSBkbyBwcm9jZXNzbyBzZW0NCnBlcmRlciBwcmVjaXPDo28uIEFtYm9zIGFwcmVzZW50YW0sIGNvbnR1ZG8gdW0gUk1TRSBzZW1lbGhhbnRlIGVtIHJlbGHDp8OjbyBhbw0KcHJvY2Vzc28gY29udMOtbnVvOg0KDQpgYGB7ciByZXN1bHRzPSdob2xkJ30NClJNU0UgPSBmdW5jdGlvbih4LHkpIHNxcnQobWVhbigoeC15KV4yKSkNCmNhdCgiUk1TRTpcbiIpDQpjYmluZCgiVGF1Y2hlbiIgPSBSTVNFKHpjLHp0KSwgIlJvdXdlbmhvcnN0IiA9IFJNU0UoemMsenIpKQ0KYGBgDQoNCiMjIyAgUXVlc3TDo28gNA0KDQojIyMjIEVzdGltZSBwcm9jZXNzb3MgQVIoMSkgY29tIGJhc2Ugbm9zIGRhZG9zIHNpbXVsYWRvcywgdGFudG8gYSBwYXJ0aXIgZG8gVGF1Y2hlbiBxdWFudG8gbyBkZSBSb3V3ZW5ob3JzdC4gUXXDo28gcHLDs3hpbW8gZWxlcyBlc3TDo28gZG8gcHJvY2Vzc28gZ2VyYWRvciBkZSBkYWRvcyByZWFsPyBTZSBlbGVzIG7Do28gZXN0aXZlcmVtIG11aXRvIHByw7N4aW1vcywgdXRpbGl6ZSBtYWlzIHBvbnRvcy4NCg0KDQoNCmBgYHtyfQ0KZml0X3p0ID0gbG0oenQgfiBzaGlmdCh6dCkgKyAwKQ0KZml0X3pyID0gbG0oenIgfiBzaGlmdCh6cikgKyAwKQ0KDQpjb2VmcyA9IGMoZml0X3p0JGNvZWZmaWNpZW50cywNCiAgICAgICAgICBmaXRfenIkY29lZmZpY2llbnRzKQ0KDQpzaWdtYXMgPSBjKHN1bW1hcnkoZml0X3p0KSRzaWdtYSwNCiAgICAgICAgICAgc3VtbWFyeShmaXRfenIpJHNpZ21hKQ0KDQp0YmwgPSBkYXRhLmZyYW1lKG1ldGhvZCA9IGMoIlRhdWNoZW4gOSBwdHMiLCANCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAiUm91d2VuaG9yc3QgOSBwdHMiKSwgDQogICAgICAgICAgICAgICAgIHJobyA9IGNvZWZzLCBzaWdtYSA9IHNpZ21hcykNCnRibA0KYGBgDQoNClBlbG8gbcOpdG9kbyBkZSBUYXVjaGVuIGNvbSA5IHBvbnRvcyBubyBncmlkLCBlc3RpbWFtb3MgJFxoYXR7XHNpZ21hfSA9IDAuMDA4JCwNCnZhbG9yIHJlbGF0aXZhbWVudGUgZGlzdGFudGUgZG8gJFxzaWdtYSQgdmVyZGFkZWlyby4gQXVtZW50YW5kbyBvIG7Dum1lcm8gZGUgDQpwb250b3Mgbm8gZ3JpZCBwYXJhIDUxLCBvYnRlbW9zIHVtYSBlc3RpbWF0aXZhIHNlbWVsaGFudGUgw6Agb2J0aWRhIHBhcmEgbw0KbcOpdG9kbyBkZSBSb3V3ZW5ob3JzdCBjb20gOSBwb250b3MuDQoNCmBgYHtyfQ0KdGhldGFUID0gdGF1Y2hlbl9ncmlkKE49NTEsIG09Mywgc2lnbWEsIHJobykNClBUID0gdGF1Y2hlbl9QKHRoZXRhVCwgcmhvLCBzaWdtYSwgbXUpDQp6dGwgPSBnZW5fcHJvY2Vzcyh0aGV0YVQsIFBULCBlcHMsIG11LCBzaWdtYSkNCg0KZml0X3p0bCA9IGxtKHp0bCB+IHNoaWZ0KHp0bCkgKyAwKQ0KcmJpbmQodGJsLCBjYmluZChtZXRob2QgPSAiVGF1Y2hlbiAyNSBwdHMiLCByaG8gPSBmaXRfenRsJGNvZWZmaWNpZW50cywgc2lnbWEgPSBzdW1tYXJ5KGZpdF96dGwpJHNpZ21hKSkNCmBgYA0KDQoNCg0KDQojIyMgIFF1ZXN0w6NvIDUNCg0KVmFtb3MgcmVmYXplciBvcyBleGVyY8OtY2lvcyBhY2ltYSBjb20gJFxyaG8gPSAwLjk5JC4NCg0KYGBge3IgZWNobz1GQUxTRX0NCnJobyA9IDAuOTkNCk4gPSA5DQoNCiMgUHJvY2Vzc28gY29udMOtbnVvDQp6YyA9IHJlcCgwLCB0aW1lX3NwYW4pDQpmb3IodCBpbiAyOnRpbWVfc3BhbikNCiAgemNbdF0gPSByaG8gKiB6Y1t0LTFdICsgZXBzW3RdDQoNCiMgVGF1Y2hlbg0KdGhldGFUID0gdGF1Y2hlbl9ncmlkKE4sIG0sIHNpZ21hLCByaG8pDQpQVCA9IHRhdWNoZW5fUCh0aGV0YVQsIHJobywgc2lnbWEsIG11KQ0KenQgPSBnZW5fcHJvY2Vzcyh0aGV0YVQsIFBULCBlcHMsIG11LCBzaWdtYSkNCg0KdGJsID0gZGF0YS5mcmFtZSh0ID0gMTp0aW1lX3NwYW4sIHpjLCB6dCkNCg0KcDEgPSBnZ3Bsb3QodGJsLCBhZXModGJsJHQpKSArDQogIGdlb21fbGluZShhZXMoeSA9IHRibCR6YywgY29sID0gIkNvbnRpbnVvdXMiKSwgc2l6ZSA9IDAuMikgKw0KICBnZW9tX2xpbmUoYWVzKHkgPSB0YmwkenQsICBjb2wgPSAiVGF1Y2hlbiBEaXNjcmV0aXphdGlvbiIpLCBzaXplID0gMC41KSArDQogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IHRoZXRhVCwgc2l6ZSA9IDAuMSwgbGluZXR5cGUgPSAzLCBhbHBoYSA9IDAuNSkgKw0KICBzY2FsZV95X2NvbnRpbnVvdXMoYnJlYWtzID0gcm91bmQodGhldGFULCA0KSwgbGltaXRzID0gYygtLjE1LCAuMTUpKSArDQogIHNjYWxlX2NvbG9yX21hbnVhbChuYW1lID0gIlByb2Nlc3MiLCB2YWx1ZXMgPSBjKCJDb250aW51b3VzIj0gMSwgIlRhdWNoZW4gRGlzY3JldGl6YXRpb24iPSA0KSkgKw0KICB0aGVtZV9jbGFzc2ljKCkgKw0KICB5bGFiKCJ6IikgKw0KICB4bGFiKCJ0IikNCg0KdGJsMiA9IHRibFsxOjUwMCwgXQ0KDQpwMiA9IGdncGxvdCh0YmwyLCBhZXModGJsMiR0KSkgKw0KICBnZW9tX2xpbmUoYWVzKHkgPSB0YmwyJHpjLCBjb2wgPSAiQ29udGludW91cyIpLCBzaXplID0gMC4yKSArDQogIGdlb21fbGluZShhZXMoeSA9IHRibDIkenQsICBjb2wgPSAiVGF1Y2hlbiBEaXNjcmV0aXphdGlvbiIpLCBzaXplID0gMC41KSArDQogIGdlb21faGxpbmUoeWludGVyY2VwdCA9IHRoZXRhVCwgc2l6ZSA9IDAuMSwgbGluZXR5cGUgPSAzLCBhbHBoYSA9IDAuNSkgKw0KICBzY2FsZV95X2NvbnRpbnVvdXMoYnJlYWtzID0gcm91bmQodGhldGFULCA0KSwgbGltaXRzID0gYygtLjE1LCAuMTUpKSArDQogIHNjYWxlX2NvbG9yX21hbnVhbChuYW1lID0gIlByb2Nlc3MiLCB2YWx1ZXMgPSBjKCJDb250aW51b3VzIj0gMSwgIlRhdWNoZW4gRGlzY3JldGl6YXRpb24iPSA0KSkgKw0KICB0aGVtZV9jbGFzc2ljKCkgKw0KICB5bGFiKCJ6IikgKw0KICB4bGFiKCJ0IikNCg0KcGwxID0gZ2dwbG90bHkocDEsIHRvb2x0aXAgPSBOVUxMKQ0KcGwyID0gZ2dwbG90bHkocDIsIHRvb2x0aXAgPSBOVUxMKQ0KDQpwbCA9IHN1YnBsb3QocGwxLCBwbDIsIG5yb3dzID0gMikNCnBsJGhlaWdodCA9IDYwMA0KcGwkd2lkdGggPSAxLjUgKiBwbCRoZWlnaHQNCg0KcGwNCmBgYA0KDQpgYGB7ciBlY2hvPUZBTFNFfQ0KTiA9IDkNCnRoZXRhUiA9IHRhdWNoZW5fZ3JpZChOLCBzcXJ0KE4tMSksIHNpZ21hLCByaG8pDQp6ciA9IGdlbl9wcm9jZXNzKHRoZXRhUiwgcm91d2VuX1AoTiwgcmhvKSwgZXBzLCBtdSwgc2lnbWEpDQoNCnRibCA9IGRhdGEuZnJhbWUodCA9IDE6dGltZV9zcGFuLCB6YywgenIpDQoNCnAxID0gZ2dwbG90KHRibCwgYWVzKHRibCR0KSkgKw0KICBnZW9tX2xpbmUoYWVzKHkgPSB0YmwkemMsIGNvbCA9ICJDb250aW51b3VzIiksIHNpemUgPSAwLjIpICsNCiAgZ2VvbV9saW5lKGFlcyh5ID0gdGJsJHpyLCAgY29sID0gIlRhdWNoZW4gRGlzY3JldGl6YXRpb24iKSwgc2l6ZSA9IDAuNSkgKw0KICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSB0aGV0YVIsIHNpemUgPSAwLjEsIGxpbmV0eXBlID0gMywgYWxwaGEgPSAwLjUpICsNCiAgc2NhbGVfeV9jb250aW51b3VzKGJyZWFrcyA9IHJvdW5kKHRoZXRhUiwgNCksIGxpbWl0cyA9IGMoLS4xNSwgLjE1KSkgKw0KICBzY2FsZV9jb2xvcl9tYW51YWwobmFtZSA9ICJQcm9jZXNzIiwgdmFsdWVzID0gYygiQ29udGludW91cyI9IDEsICJUYXVjaGVuIERpc2NyZXRpemF0aW9uIj0gImJyb3duMyIpKSArDQogIHRoZW1lX2NsYXNzaWMoKSArDQogIHlsYWIoInoiKSArDQogIHhsYWIoInQiKQ0KDQp0YmwyID0gdGJsWzE6NTAwLCBdDQoNCnAyID0gZ2dwbG90KHRibDIsIGFlcyh0YmwyJHQpKSArDQogIGdlb21fbGluZShhZXMoeSA9IHRibDIkemMsIGNvbCA9ICJDb250aW51b3VzIiksIHNpemUgPSAwLjIpICsNCiAgZ2VvbV9saW5lKGFlcyh5ID0gdGJsMiR6ciwgIGNvbCA9ICJUYXVjaGVuIERpc2NyZXRpemF0aW9uIiksIHNpemUgPSAwLjUpICsNCiAgZ2VvbV9obGluZSh5aW50ZXJjZXB0ID0gdGhldGFSLCBzaXplID0gMC4xLCBsaW5ldHlwZSA9IDMsIGFscGhhID0gMC41KSArDQogIHNjYWxlX3lfY29udGludW91cyhicmVha3MgPSByb3VuZCh0aGV0YVIsIDQpLCBsaW1pdHMgPSBjKC0uMTUsIC4xNSkpICsNCiAgc2NhbGVfY29sb3JfbWFudWFsKG5hbWUgPSAiUHJvY2VzcyIsIHZhbHVlcyA9IGMoIkNvbnRpbnVvdXMiPSAxLCAiVGF1Y2hlbiBEaXNjcmV0aXphdGlvbiI9ICJicm93bjMiKSkgKw0KICB0aGVtZV9jbGFzc2ljKCkgKw0KICB5bGFiKCJ6IikgKw0KICB4bGFiKCJ0IikNCg0KcGwxID0gZ2dwbG90bHkocDEsIHRvb2x0aXAgPSBOVUxMKQ0KcGwyID0gZ2dwbG90bHkocDIsIHRvb2x0aXAgPSBOVUxMKQ0KDQpwbCA9IHN1YnBsb3QocGwxLCBwbDIsIG5yb3dzID0gMikNCnBsJGhlaWdodCA9IDYwMA0KcGwkd2lkdGggPSAxLjUgKiBwbCRoZWlnaHQNCg0KcGwNCmBgYA0KDQoNCg0KTyBtw6l0b2RvIGRlIFJvdXdlbmhvcnN0IGZ1bmNpb25hIG1lbGhvciBwYXJhICRccmhvID0gMC45OSQgY29tIDkgcG9udG9zIG5vIGdyaWQuDQpPYnRlbW9zIG9zIHNlZ3VpbnRlcyB2YWxvcmVzIGRlICRcaGF0e1xyaG99JCBlICRcaGF0e1xzaWdtYX0kOg0KDQoNCmBgYHtyIGVjaG89RkFMU0V9DQpmaXRfenQgPSBsbSh6dCB+IHNoaWZ0KHp0KSArIDApDQpmaXRfenIgPSBsbSh6ciB+IHNoaWZ0KHpyKSArIDApDQoNCmRhdGEuZnJhbWUobWV0aG9kID0gYygiVGF1Y2hlbiA5IHB0cyIsICJSb3V3ZW5ob3JzdCA5IHB0cyIpLA0KICAgICAgInJobyIgPSBjKGZpdF96dCRjb2VmZmljaWVudHMsIGZpdF96ciRjb2VmZmljaWVudHMpLA0KICAgICAgInNpZ21hIiA9IGMoc3VtbWFyeShmaXRfenQpJHNpZ21hLCBzdW1tYXJ5KGZpdF96cikkc2lnbWEpKQ0KYGBgDQoNCg==